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ABSTRACT 

The  problem  of  localising  a  point  source  of  gamma  radiation  is  considered.  A 
simplified  analytical  approach  based  on  the  inverse  distance  square  law  as  well 
as  several  probabilistic  approaches  are  described.  The  problem  is  studied  using 
the  Cramer-Rao  bound  (CRB)  analysis,  which  quantifies  the  accuracy  with 
which  it  is  possible  to  localise  the  source  and  estimate  its  activity.  Simulated 
and  real  radiological  survey  data  are  used  to  investigate  the  performance  of 
the  algorithms. 
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Radiological  Source  Localisation 

Executive  Summary 

The  increasing  threat  of  chemical,  biological  and  radiological  (CBR)  attacks  has  re¬ 
sulted  in  a  significant  interest  in  research  on  countering  such  attacks.  Our  research  focuses 
on  countering  radiological  attacks,  which  may,  for  example,  be  carried  out  using  radiolog¬ 
ical  dispersion  devices  or  dirty  bombs.  The  ability  to  rapidly  localise  a  radiological  source 
can  assist  emergency  responders  to  disable,  isolate  or  safely  remove  such  a  device. 

This  report  describes  some  preliminary  work  we  have  carried  out  in  the  area  of  radi¬ 
ological  source  modelling  and  localisation.  This  work  concerned  localisation  of  a  single 
fixed  gamma  radiation  source  of  unknown  activity  level.  The  accuracy  with  which  the 
source  location  and  the  activity  could  be  estimated  was  studied  using  the  Cramer  Rao 
bound  analysis.  A  simple  deterministic  analytical  approach  as  well  as  several  probabilistic 
estimation  techniques  were  investigated  using  simulated  and  real  measurement  data. 

The  inverse  square  law-based  deterministic  solution  developed  in  this  work  uses  radi¬ 
ation  measurements  collected  at  four  arbitrary  points  to  estimate  the  source  position  and 
activity.  This  algorithm  was  able  to  provide  reasonable  source  estimates  based  on  real 
data  collected  using  the  Low  Cost  Advanced  Airborne  Radiological  Survey  (LCAARS) 
system  developed  by  DSTO. 

The  maximum  likelihood  estimator  and  a  nonlinear  least  squares  estimator  yielded 
quite  accurate  estimates.  While  the  maximum  likelihood  is  an  asymptotically  efficient 
estimator,  it  is  a  batch  algorithm  and  is  unattractive  for  operational  use.  An  inexact 
recursive  least  squares  algorithm  was  developed  and  it  produced  good  estimates  when 
applied  to  real  data.  Unscented  Kalman  filter  (UKF)  and  extended  Kalman  filter  (EKF) 
algorithms  were  also  investigated.  The  UKF  approach  performed  well  but  the  FKF  was 
divergent  and  could  not  return  acceptable  source  estimates. 

The  preliminary  investigations  described  in  this  report  have  provided  useful  insights 
to  the  issues  associated  with  radiological  source  estimation.  Future  work  will  investigate 
more  robust  estimation  techniques  and  also  consider  multiple  possibly  moving  sources, 
localisation  in  the  presence  of  obstructions  and  observer  motion  optimisation. 


DSTO-TR-1988 


DSTO-TR-1988 


Authors 


Ajith  Gunatilaka 

Human  Performance  and  Protection  Division 

Ajith  Gunatilaka  obtained  his  PhD  in  Electrical  Engineering 
from  The  Ohio  State  University,  USA  in  2000.  He  obtained  a 
B.Sc.  (Engineering)  degree  with  First  Class  Honours  in  Elec¬ 
tronics  and  Telecommunication  Engineering  from  The  Univer¬ 
sity  of  Moratuwa,  Sri  Lanka  in  1990.  Dr  Gunatilaka  joined 
DSTO’s  HPPD  in  November  2005  as  a  defence  scientist. 


Branko  Ristic 

Information  Surveillance  and  Reconnaissance  Division 

Branko  Ristic  has  been  involved  in  R&D  projects  related  to  sig¬ 
nal  processing,  estimation,  tracking  and  classification  for  more 
than  20  years.  Dr.  Ristic  received  all  of  his  degrees  in  electrical 
engineering:  Ph.D.  from  QUT  1995,  M.Sc.  Belgrade  Univer¬ 
sity  1991,  B.  Eng.  University  of  Novi  Sad  1984.  Since  1996 
he  has  been  with  DSTO,  where  his  role  has  been  to  carry  ont 
research,  develop  new  capabilities  and  provide  technical  advice 
to  the  ADO  on  the  topics  of  target  tracking,  sensor  fnsion,  and 
network  enabled  snrveillance. 

Dr  Ristic  has  pnblished  over  80  technical  papers  and  presented 
several  invited  talks  and  short  conrses  in  Australia,  Europe  and 
the  US.  He  is  lecturing  a  post-graduate  subject  at  Adelaide 
University  (snbject  ’’Mnlti-Sensor  Data  Fusion”).  He  served 
on  technical  committees  of  several  international  conferences, 
and  is  the  Ghair  of  the  Fonrth  Australian  Data  Fusion  Sym- 
posinm  (IDG-07).  Dr  Ristic  won  twice  the  best  paper  award 
(at  Information  Fusion  2005,  Philadelphia,  USA  and  DIGTA 
2005,  Gairns,  Australia),  and  recently  co-anthored  a  best-selling 
book:  Beyond  the  Kalman  filter:  Particle  filters  for  tracking  ap¬ 
plications  (Artech  House,  2004). 


v 


DSTO-TR~1988 


Ralph  Gailis 

Human  Performance  and  Protection  Division 

Ralph  Gailis  completed  his  undergraduate  studies  in  physics 
and  mathematics  (Honours)  at  the  University  of  Melbourne  in 
1992.  He  went  on  to  do  a  Ph.D.  on  the  topic  of  “Plasma  Physics 
in  the  Early  Universe”  in  the  School  of  Physics,  University  of 
Melbourne  from  1993-1996.  1997-1998  saw  him  continue  on 
in  the  same  department  as  a  postdoctoral  research  fellow,  un¬ 
til  his  recruitment  in  the  Combatant  Protection  and  Nutrition 
Branch  at  DSTO  in  1999.  Research  interests  at  the  University 
of  Melbourne  included  statistical  mechanics,  general  relativity, 
cosmology  and  plasma  physics. 

Since  working  at  DSTO,  Ralph  has  developed  a  programme 
on  hazard  assessment  of  chemical,  biological  and  radiological 
weapons.  The  main  focus  of  scientific  research  in  this  area 
has  been  in  atmospheric  dispersion  modelling  and  turbulence 
theory.  This  has  involved  some  challenging  theoretical  work,  as 
well  as  experimental  science  in  wind  tunnel  and  water  channel 
simulations  of  dispersion  processes.  The  work  has  also  involved 
operational  support  to  provide  the  ADF  with  current  state-of- 
the-art  CBRN  hazard  prediction  models,  which  have  been  used 
for  a  number  of  high  profile  events. 

Ralph  spent  over  a  year  in  2002-2003  working  at  DRDC  (Suffield), 
Alberta,  Canada  under  a  Defence  Science  Fellowship,  greatly 
extending  his  knowledge  on  the  science  behind  CBR  hazard 
assessment  and  dispersion  modelling.  He  has  also  served  the 
last  six  years  as  National  Leader,  TTCP  CBD  Group  Technical 
Panel  9  (Hazard  Assessment).  These  invaluable  international 
interactions  have  helped  the  development  of  the  DSTO  pro¬ 
gramme  in  this  area  immensely. 


VI 


DSTO-TR~1988 


Contents 


1  Introduction  1 

2  Description  of  the  problem  1 

3  Deterministic  solutions  3 

4  Statistical  data  modelling  5 

4.1  Dose  rate  data .  6 

4.2  Count  rate  data  .  9 

5  Probabilistic  solutions  11 

6  Cramer-Rao  analysis  13 

6.1  Derivation .  13 

6.2  Analysis .  14 

7  Estimation  algorithms  and  their  performance  14 

7.1  Algorithms .  14 

7.1.1  Maximum  likelihood  estimation .  14 

7.1.2  Least  squares  (LS)  with  reduced  search  dimension .  16 

7.1.3  An  approximate  recursive  LS  algorithm .  17 

7.1.4  EKF  and  UKF  .  18 

7.1.5  EKF  equations  .  18 

7.1.6  UKF  equations .  19 

7.2  Simulation  results  .  20 

8  Application  to  real  radiological  survey  data  20 

8.1  Holsworthy  field  trial  data .  20 

8.2  Deterministic  solutions .  22 

8.3  EKF,  UKF  and  MLE  .  24 

8.4  LS  with  reduced  search  dimension .  25 

8.5  The  approximate  recursive  LS  algorithm .  31 


vii 


9  Conclusions 


32 


DSTO-TR-1988 


References 


36 


Appendices 


A  Derivation  of  the  analytical  solution 


37 


viii 


DSTO-TR~1988 


Figures 

1  Geometry  of  the  radiological  point  source  localisation  problem  .  2 

2  Raw  dose  rate  data  measured  with  the  gamma-beta  probe .  6 

3  Discreteness  of  dose  rate  data  measured  with  the  gamma-beta  probe .  7 

4  The  histogram  for  normalised  dose  rate  data  and  the  Poisson  and  Gaussian 

distribution  fitting .  7 

5  Data  fitting  to  dose  rate  data  measured  with  the  gamma-beta  probe .  9 

6  Raw  count  rate  data  collected  with  the  Micro-R  probe .  10 

7  Normalised  histograms  and  the  Gaussian  fit  for  Micro-R  probe  data .  10 

8  The  characterisation  of  Micro-R  probe  data .  11 

9  GRB  analysis:  the  error  standard  deviation  in  (a)  position  and  (b)  activity  .  15 

10  The  scenario  used  for  simulation  .  21 

11  The  estimation  error  in  position .  21 

12  The  estimation  error  in  source  activity .  21 

13  An  example  application  of  the  four  point  based  deterministic  source  estima¬ 
tion  algorithm .  23 

14  The  statistical  spread  of  the  estimated  source  location  over  a  set  of  30  four- 

point  inputs . 23 

15  The  normalised  histogram  for  background  dose  rate  data .  24 

16  Trial  data  analysis .  25 

17  Measured  dose  rates  and  Afc(x)  corresponding  to  the  MLE .  26 

18  The  squared  error  surface  obtained  for  ^^^Gs  count  rate  data . 27 

19  Gomparison  of  ^^^Gs  count  rate  measurements  and  prediction  based  on  LS 

estimation .  27 

20  The  squared  error  surface  obtained  for  ^^^Gs  dose  rate  data . 28 

21  Gomparison  of  ^^'^Gs  dose  rate  measurements  and  prediction  based  on  LS 

estimation .  28 

22  The  squared  error  surface  obtained  for  ®*^Go  count  rate  data . 29 

23  Gomparison  of  ®^Go  count  rate  measurements  and  prediction  based  on  LS 

estimation .  30 

24  The  squared  error  surface  obtained  for  ®*^Go  dose  rate  data  . 30 

25  Gomparison  of  ®*^Go  dose  rate  measurements  and  prediction  based  on  LS 

estimation .  31 

26  The  squared  error  surface  obtained  for  background  measurements . 32 

27  The  comparison  of  the  batch  LS  algorithm  with  the  approximate  recursive 

LS  algorithm  for  the  case  of  ^^^Gs  count  rate  data  .  33 


IX 


DSTO-TR-1988 


Tables 


1  The  mean  value  and  the  standard  deviation  of  dose  rate  data .  8 

2  Mean  value  and  standard  deviation  of  count  rate  data .  9 

3  An  example  set  of  four  measurements  from  the  ^^^Cs  count  rate  data  used  to 

demonstrate  source  estimation  using  the  four  point  estimation  algorithm  ...  22 


X 


DSTO-TR~1988 


1  Introduction 


Enhancing  the  chemical,  biological  and  radiological  (CBR)  situational  awareness  of  the 
battlefield  is  an  important  goal  of  DSTO  research.  Accurate  assessment  and  modelling  of 
CBR  hazards  are  essential  for  the  purposes  of  safely  locating  and  removing  or  isolating 
them.  Hazard  modelling  also  makes  it  possible  to  predict  the  impact  of  CBR  hazards  on 
surrounding  geographical  areas  and  populations.  Atmospheric  dispersion  modelling  tools, 
such  as  Hazard  Prediction  and  Assessment  Capability  (HPAC),  together  with  accurate 
weather  data  provide  a  powerful  and  convenient  way  to  determine  how  a  release  of  a 
known  CBR  source  would  disperse  and  what  kind  of  hazardous  exposure  it  would  cause  in 
an  area.  This  forward  modelling  capability  provided  by  HPAC  and  similar  tools  is  useful 
in  assessing  who  is  vulnerable  to  the  effects  of  a  known  CBR  source  that  exists  at  a  known 
location.  It  would  generally  be  expected  that  a  CBR  source  is  not  known.  Sometimes 
approximate  prior  knowledge  about  the  CBR  source  and  its  location  may  be  available,  for 
example  through  intelligence  information.  The  exact  nature  of  the  source  and  its  location 
need  to  be  ascertained  by  using  sensor  measurements  and  appropriate  source  estimation 
algorithms.  A  previous  DSTO  report[l]  outlined  DSTO’s  planned  work  towards  CBR 
source  term  estimation,  hazard  modelling  and  data  fusion. 

Because  CB  detection  technologies  are  relatively  immature  and  also  because  other 
nations  in  the  TTCP  have  devoted  limited  attention  to  the  radiological  problem  [1] ,  it  has 
been  decided  to  focus  DSTO’s  initial  work  in  this  area  to  radiological  source  estimation. 
Another  reason  for  choosing  to  initially  work  with  the  radiological  source  problem  is  the 
availability  of  the  DSTO-developed  Low  Cost  Advanced  Radiological  Survey  (LCAARS) 
system  for  data  collection  to  support  testing  and  validation  of  algorithms.  DSTO  intends 
to  undertake  work  in  CB  hazard  modelling  and  data  fusion  in  future.  This  report  describes 
our  initial  work  in  radiological  hazard  modelling  and  source  term  estimation. 

The  ultimate  goal  of  radiological  source  estimation  research  is  to  be  able  to  track  a 
dispersing  cloud  of  radioactive  particles  and  predict  the  integrated  radiation  dose  emitted 
in  all  directions  from  the  cloud.  Before  we  can  solve  this  rather  complex  problem,  it 
is  essential  to  gain  a  better  understanding  of  the  relevant  issues  by  first  tackling  the 
relatively  simpler  problem  of  estimating  a  single  radiological  point  source.  Once  this 
initial  problem  is  tackled,  future  research  will  consider  the  problem  of  estimating  multiple 
point  gamma  sources.  The  overview  of  research  plan  discussed  in  Section  4.2  of  [1]  suggests 
to  consider  the  problem  of  locating  a  point  gamma  radiation  source  as  a  first  enhancement 
of  LCAARS,  and  to  test  these  algorithms  in  a  field  trial. 


2  Description  of  the  problem 

Radiation  is  energy  that  originates  in  atomic  or  nuclear  processes  and  travels  through 
space  and  materials.  Radiation  may  include  fast  electrons  including  beta  particles,  heavy 
charged  particles  such  as  alpha  particles,  electromagnetic  radiation  including  X-rays  and 
gamma  rays,  and  neutron  radiation  [2].  Alpha  radiation  travels  only  short  distances 
through  air  and  is  unable  to  penetrate  clothing  or  human  skin.  Alpha-emitting  material 
may,  however,  be  hazardous  if  they  enter  the  body  through  inhalation  or  other  means. 
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Beta  radiation  has  higher  penetration  than  alpha  radiation  and  may  moderately  penetrate 
human  skin.  Internal  deposition  of  beta-emitting  material  can  be  harmful.  Devices  such 
as  thin-window  Geiger-Miiller  probes  can  detect  the  presence  of  alpha  or  beta  sources. 
Because  of  limited  distances  travelled  by  alpha  and  beta  radiation,  stand-off  detection  of 
these  sources  is  difficult.  Gamma  rays,  in  contrast,  are  highly  penetrating  and  may  be 
detected  at  large  distances  from  the  source. 

Developing  a  capability  to  accurately  pinpoint  the  location  of  an  unknown  source 
of  gamma  radiation,  (eg.  a  dirty  bomb)  can  help  rapid  identification  and  removal  of 
such  threats.  As  mentioned  above,  currently,  radiation  survey  data  collected  with  the 
LGAARS  system  has  to  be  loaded  into  a  mapping  software  and  a  two-dimensional  intensity 
plot  produced  to  visually  identify  the  area  with  the  highest  intensity  as  the  putative 
source  location.  If  the  source  lies  outside  the  surveyed  area,  then  its  position  cannot  be 
determined  using  this  approach.  The  gradient  of  the  data  may  provide  an  indication  of 
the  direction  of  the  source  location,  which  may  then  be  used  as  a  cue  for  a  fresh  survey. 
The  goal  of  the  current  work  is  to  investigate  approaches  to  improve  the  accuracy  and  the 
efficiency  of  locating  gamma  radiation  sources,  preferably  in  real  time. 

Gonsider  a  single  radiological  point  source  located  at  {xo,yo)  in  the  two-dimensional 
Gartesian  coordinate  system  (Fig.  1).  For  simplicity,  the  elevation  of  the  source  is  not 
considered.  A  radiation  survey  instrument  is  used  to  collect  measurements  of  radiation 
Zk  at  Gartesian  coordinates  {xk,yk),  where  k  =  1, . . . ,  and  N  is  the  total  number  of 
measurements.  It  is  assumed  that  {{xk,  yk)',k  =  1, . . . ,  N}  are  known  (measured  using  a 
GPS  receiver  for  example).  The  goal  is  to  estimate  the  location  and  the  strength  of  the 
unknown  source  based  on  these  measurements. 


Figure  1:  Geometry  of  the  radiologieal  point  souree  loealisation  problem 


The  dose  rate  D  (measured  in  units  of  gSv/h)  at  a  distance  r  (m)  from  a  point  gamma 
source  can  be  written  as  [3]: 

^  ME 

-  (1) 

where  M  is  the  source  activity  in  MBq  and  E  is  the  gamma  energy  per  disintegration  in 
MeV. 

We  can  rewrite  Equation  (1)  in  the  following  familiar  form  of  the  inverse  squared 
distance  law: 
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where 


ME 

~Q~ 


(2) 

(3) 


Here  the  single  parameter  I  is  used  to  characterise  the  product  of  the  radioactivity 
of  the  source  and  the  gamma  energy  per  disintegration.  From  Equation  (2),  I  can  also 
be  interpreted  as  the  dose  rate  measured  at  a  distance  of  1  m  from  the  source.  In  some 
cases  the  gamma  energy  E  may  be  known,  for  example  through  gamma  spectroscopy  and, 
therefore,  an  estimate  of  I  can  be  used  to  compute  the  source  activity  M.  Therefore,  for 
simplicity,  in  this  report  we  may  refer  to  I  as  activity. 

When  the  LCAARS  system  is  used  in  conjunction  with  the  standard  gamma  probe,  it 
measures  the  dose  rate  and  records  in  units  of  mSv/h.  Therefore  the  estimated  activity 
will  be  in  GBq. 

If  the  LCAARS  system  is  used  with  the  Micro-R  probe,  the  survey  data  will  be  count 
rates  expressed  as  counts  per  minute  (CPM).  In  this  case,  Equation(l)  must  be  modified 
to  incorporate  the  energy  response  (CPM//iSv/h)  of  the  probe.  Therefore, 


Count  Rate 


ME 

/ 


X  Energy  Response 


(4) 

(5) 


where 

ME 

I  =  -  X  Energy  Response  (6) 

6 

From  Equation  (5)  we  can  interpret  I  as  the  count  rate  measured  at  a  distance  of  1  m 
away  from  the  source. 


The  Micro-R  probe  has  an  energy  response  of  12200  CPM/^Sv/h  for  0.662  MeV  gamma 
radiation  of  ^^'^Cs  and  6700  CPM/^uSv/h  for  1.25  MeV  average  gamma  energy  of  ®^Co 
source  [4].  Assuming  that  the  energy  of  the  measured  radiation  is  known,  the  source 
radiative  activity  M  can  be  computed  using  Equation  (6)  and  the  estimated  value  of  /. 
Therefore,  we  can  represent  the  unknown  point  gamma  radiation  source  by  defining  a 
parameter  vector  as  follows: 

x=[xo  2/0  lY  (7) 

where  T  denotes  the  matrix  transpose. 


Our  goal  is  to  estimate  vector  x  using  measurements  {zk]k  =  1,...,^}  collected  at 
coordinates  {(x^,  2/fc);  A:  =  1, . . . ,  N}. 


3  Deterministic  solutions 


In  simple  terms,  the  point  source  estimation  problem  can  be  considered  as  estimating 
the  three  unknown  parameters  {xo,yQ,I).  In  principle,  only  three  equations  in  terms  of 
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these  three  unknowns  are  necessary  to  obtain  an  exact  solution.  Because  the  dose  rate 
or  the  count  rate  Zk  measured  at  any  location  {xk,yk)  depends  on  the  three  unknown  pa¬ 
rameters  through  the  inverse  square  law,  measurements  obtained  at  three  different  points 
are  all  that  is  required  to  estimate  the  source  location  and  strength.  Because  the  solution 
of  these  equations  involves  several  square  root  operations,  it  produces  several  spurious 
solutions  in  addition  to  the  true  estimate,  making  this  solution  difficult  to  use.  By  includ¬ 
ing  the  measurement  from  a  fourth  point,  these  square  root  operations  and,  hence,  the 
spurious  solutions  can  be  avoided. 

Let  us  consider  measurements  {z^, /c  =  1, 2, 3, 4}  collected  at  four  arbitrary  points 
{xk-iVk)-  If  we  disregard  noise  and  background  radiation,  these  measurements  can  be 
expressed  as: 

Zk{y^)=I/rl,  A:  =  1,2, 3, 4  (8) 


where 

rk  =  \J {xk  -  xo)^  +  (2/fc  -  yoP  (9) 

is  the  distance  from  the  source  (xo,yo)  to  the  kth.  measurement  point  {xk,yk)- 

As  shown  in  Appendix  A,  Equation  (8)can  be  solved  analytically  for  xq,  yo  and  /, 
which  gives: 


(10) 

(11) 

(12) 


where 

Ai  = 

K2  = 

Li  =  {xl  -  xl)  +  {yf  -  yl) 

L2  =  {xl  —  xl)  +  {yl  —  yl) 

As  =  {xl  -  xl)  -f  {yl  -  yl) 

ai  =  2Ki{xi  -  x^)  -  2{xi  -  X2) 

bi  =  2Ki{yi  -  yi)  -  2{yi  -  y2) 
Cl  =  Li  —  K1L2 

02  =  2K2{xi  -  X4,)  -  2{xi  -  X2) 

62  =  2K2{yi  -  y^)  -  2(yi  -  2/2) 

C2  =  Li  —  K2L3. 


Xq  = 

yo  = 

I  = 


biC2  -  Cl 62 
0162  -  6102 

aiC2  —  Cl  02 


and 


61 02  -  0162 

zi  [(xi  -  xo)^  -f  {yi  -  yo)^] 
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Since  the  radioactive  decay  is  a  stochastic  process,  the  inverse  distance  square  law  is 
applicable  to  the  measurements  only  in  an  average  sense.  Therefore,  the  count  rate  or 
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the  dose  rate  measurements  should  truly  be  the  averages  taken  over  a  sufficiently  large 
number  of  measurements  collected  at  each  measurement  point. 

The  above  analytical  solution  based  on  four  measurements  is  quite  simple  and  may  be 
useful  as  a  tool  to  rapidly  obtain  at  least  a  crude  estimate  of  a  radiological  point  source. 
There  are,  however,  many  limitations  in  the  applicability.  The  solution  assumes  a  single 
point  source  in  open  space  and,  therefore,  not  applicable  to  multiple  source  situations  or 
when  there  are  obstacles  between  the  source  and  the  measurement  points.  For  simplicity, 
the  background  radiation  level  was  ignored  in  the  derivation  discussed  above.  This  is  not 
a  major  limitation  because  it  is  possible  in  real  measurements  to  subtract  the  average 
background  from  all  measurements,  provided  the  background  level  is  relatively  constant. 
The  results  obtained  by  applying  the  four  point  estimation  algorithm  to  real  measurement 
data  is  discussed  in  Subsection  8.2. 


4  Statistical  data  modelling 

Accurate  statistical  modelling  of  data  is  essential  to  the  development  of  effective  es¬ 
timation  algorithms.  A  set  of  controlled  experiments  was  carried  out  in  June  2006  to 
collect  data  to  characterise  the  radiation  survey  measurements.  These  experiments  were 
conducted  in  the  radiation  laboratory  at  the  DSTO  Fishermans  Bend  site  using  the  low 
cost  advanced  airborne  radiological  survey  (LCAARS)  system. 

The  CBRN  Defence  Centre  at  DSTO  prototyped  this  survey  system  using  an  AN /PDR- 
77  radiation  survey  meter  equipped  with  an  RS232  interface  module,  a  gamma  radiation 
detector,  a  GPS  receiver,  and  survey  software.  Although  the  system  was  originally  de¬ 
signed  for  airborne  survey,  as  the  name  suggests,  the  DSTO-developed  prototype  was 
primarily  developed  for  ground  vehicle  based  survey. 

Two  different  radiation  detectors;  a  standard  gamma  probe  and  a  Micro-R  probe  were 
used.  With  the  standard  gamma  probe,  which  is  a  Geiger-Miiller  tube,  the  AN/PDR-77  is 
capable  of  measuring  gamma  radiation  over  a  wide  range  of  dose  rates  without  saturating 
at  high  levels  but  the  sensitivity  is  poor  at  low  dose  rates  close  to  background  and  at 
low  gamma  energies.  Dose  rate  data  measured  with  this  probe  were  recorded  in  milligray 
per  hour  (mGy/h).  The  Micro-R  probe  uses  a  2.54  cm  (1  inch)  diameter  and  3.81  cm 
(1.5  inch)  long  Nal  scintillation  crystal  coupled  to  a  2.54  cm  diameter  photomultiplier 
tube  operated  at  700  volts.  It  has  greater  sensitivity  to  low  energy  and  low  dose  rate 
radiation.  Measurements  with  this  probe  were  recorded  as  count  rates  in  counts  per 
minute  (CPM). 

Data  were  collected  with  each  sensor  positioned  1  m,  1.5  m,  2  m,  3  m  and  4  m  away 
from  a  1.647  GBq  ®^Co  radiation  source  which  was  enclosed  within  its  transport  container 
to  prevent  excessive  exposure  to  the  experimenter.  Background  data  in  the  absence  of  any 
radiation  sources  were  also  collected  with  each  sensor.  The  data  collection  software  was 
set  up  to  acquire  data  for  a  period  of  one  hour  for  each  different  configuration.  As  the 
system  collects  one  data  point  per  approximately  every  two  seconds,  approximately  1840 
measurements  were  collected  in  each  configuration.  The  sensors  were  stationary  during 
data  acquisition. 


5 


DSTO-TR~1988 


4.1  Dose  rate  data 

The  raw  dose  rate  data  measured  using  the  gamma  probe  (in  mSv/h)  at  different  dis¬ 
tances  from  the  source  as  well  as  the  data  collected  in  the  absence  of  a  source  are  shown 
in  Figure  2.  The  variance  of  these  data  can  be  seen  to  increase  with  the  mean  signal 
level.  Figure  3(a)  shows  the  raw  data  (zraw)  collected  3  m  away  from  the  radiation  source, 
which  show  that  the  measurements  fall  into  a  set  of  discrete  levels  including  0,  0.00021, 
0.00042,  0.00063  mSv/h  and  other  integer  multiples  of  0.00021.  A  few  measurements 
are  slightly  outside  these  levels  having  values  such  as  0.00041  or  0.00062,  which  may  be 
due  to  electronic  noise  in  the  sensor.  All  other  data  sets  also  displayed  similar  discrete¬ 
ness.  Figure  3(b)  shows  the  normalised  data  (znorm)  that  were  obtained  by  applying  the 
transformation: 


1  /  Zraw  \ 

Znorm  =  I’OUnd  -  . 

V  0.00021/ 

Figure  4  shows  the  normalised  histograms  of  each  data  set.  The  mean  value  /r  and  the 
standard  deviation  a  for  each  of  the  data  sets  are  listed  in  Table  1.  From  the  third  and 
fourth  columns  of  the  table  we  notice  that  a  ~  as  expected  for  a  Poisson  distribution. 
If  z  =  0, 1, 2, . . .  is  the  dose  rate  count,  then  the  Poisson  pmf  is  given  by  V{z]^)  = 
e~^fi^/z\.  Since  both  the  mean  and  the  variance  of  a  Poisson  distribution  are  equal 
(to  //),  for  smaller  distances  between  the  source  and  the  probe,  the  probe  readings  are 
characterised  by  higher  magnitude  but  also  by  higher  uncertainty.  A  Poisson  distribution 
with  large  jj.  is  approximately  Gaussian  with  mean  and  variance  equal  to  /x.  Observe 
from  Figure  4  that  the  Gaussian  fit  is  fairly  accurate  for  all  cases  except  for  background 
radiation. 


„  Gamma-beta  probe  data 


- Background 

Co  source  4m  away 

8  -  Co  source  3m  away 

- Co  source  2m  away 

- Co  source  1 .5m  away 


Measurement  number 


Figure  2:  Raw  dose  rate  data  measured  with  the  gamma-beta  probe 
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Gamma-beta  probe  3m  away  from  Co60  source 


(a)  Raw  data 


(b)  Normalised  data 


Figure  3:  Discreteness  of  dose  rate  data  measured  with  the  gamma-beta  probe 


Probability  distribution  of  Gamma  probe  data 


Figure  4'  The  histogram  for  normalised  dose  rate  data  and  the  Poisson  and  Gaussian 
distribution  fitting. 
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If  we  assume  that  the  mean  values  of  the  data  measured  at  different  distances  from 
the  source  to  follow  the  inverse  square  law,  then  we  can  express  these  mean  values  as  a 
function  of  distance  as  follows: 


where  is  the  value  of  fi  when  the  source  activity  /  =  0  or  the  source  distance  d  =  oo. 
Thus  Ho  represents  the  mean  value  of  the  background. 

By  plotting  the  measured  data  and  performing  a  linear  fit,  we  estimate  the  mean  value 
of  the  background  radiation  as  ~  1-5.  No  unit  is  specified  here  for  hoi  as  we  are 
considering  the  normalised  dose  rate,  which  is  dimensionless.  Although  we  assumed  the 
data  to  follow  the  inverse  square  law,  and  used  a  linear  fit  in  Figure  5(a)  accordingly,  from 
the  figure,  the  data  do  not  seem  to  follow  a  straight  line  fit  well. 

To  check  the  validity  of  this  assumption  let  us  assume  the  mean  values  of  the  data 
measured  at  different  distances  from  the  source  to  obey  the  following  rule: 


I 

^  ^  +  ^0- 

By  taking  the  logarithm  of  both  sides  of  Equation  (14)  we  obtain: 

login  -  IJ-o)  =  log  /  -  a  log  d. 


(14) 

(15) 


Figure  5(b)  shows  the  plot  of  Equation  (15),  where  no  =  0.7946  estimated  from  the 
background  data  has  been  used.  The  linear  fit  to  the  data  in  Figure  5(b)  shows  that 
a  =  — 1.8  and  log  /  =  2.9  or  /  =  exp(2.9)  18.  Therefore, 

18 

/U=  ^+0.7946.  (16) 


Because  a  =  —1.8,  rather  than  —2,  this  result  suggests  that  the  mean  values  of  the 
measured  dose  rate  data  fall  off  at  a  rate  slightly  less  than  that  predicted  by  the  inverse 
square  law.  We  believe  that  this  less-than-inverse  square  fall  off  is  due  to  scattering  from 
the  concrete  walls  of  the  laboratory.  We  assume  the  inverse  square  law  to  hold  for  radiation 
survey  data  measured  in  open  space  such  as  in  the  Holsworthy  trial.  We  intend  to  verify 
the  validity  of  this  assumption  during  the  next  field  trial  planned  for  later  this  year,  as  the 
inverse  square  law  was  a  basic  assumption  made  in  developing  the  estimation  algorithms. 


Table  1:  The  mean  value  and  the  standard  deviation  of  dose  rate  data 


Data  set 

h 

y/Jj- 

a 

1  m 

17.6339 

4.1993 

4.2565 

1.5  m 

9.7951 

3.1313 

3.1313 

2  m 

6.2687 

2.5037 

2.4900 

3  m 

3.3643 

1.8342 

1.8081 

4  m 

2.2721 

1.5074 

1.4983 

Background 

0.7946 

0.8914 

0.8726 
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Gamma  dose  rate  probe  data  Gamma  dose  rate  probe  data 


Figure  5:  Data  fitting  to  dose  rate  data  measured  with  the  gamma-beta  probe 


4.2  Count  rate  data 

A  plot  of  the  data  collected  in  the  laboratory  with  the  Micro-R  sensor  is  shown  in 
Figure  6.  The  normalised  histograms  for  these  data  are  shown  in  Figure  7.  As  in  the  case 
of  the  dose  rate  data,  here  also  we  clearly  observe  that  the  uncertainty  of  data  increases 
with  the  mean  value.  A  Gaussian  function  hts  these  data  quite  well. 

The  mean  value  fi,  ^JJi  and  the  standard  deviation  a  for  each  data  set  are  listed  in 
Table  2.  Unlike  in  the  dose  rate  data,  here  the  variance  is  not  equal  to  the  mean. 

Table  2:  Mean  value  and  standard  deviation  of  count  rate  data 


Data  set 

y/fi 

a 

1  m 

79151.303 

281.338 

3693.219 

1.5  m 

43082.790 

207.564 

2208.091 

2  m 

27020.738 

164.380 

1496.266 

3  m 

11553.149 

107.486 

865.810 

4  m 

6967.930 

83.474 

651.190 

Background 

1849.658 

43.008 

264.925 

Let  us  assume  that  cr  and  g  have  the  following  relationship: 

a  =  /3/r“. 


Therefore, 


log(cj)  =  a\og{g)  +  log(/?). 


(17) 

(18) 


By  plotting  log(cr)  versus  log(/r),  as  shown  in  Figure  8(a)  and  performing  a  linear  fit, 
we  obtain  a  =  0.69  and  j3  =  exp(0.37)  =  1.45.  Therefore,  we  can  model  the  Micro-R 
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Figure  6: 


Raw  count  rate  data  collected  with  the  Micro-R  probe 


Figure  li  Normalised  histograms  and  the  Gaussian  fit  for  Micro-R  probe  data 
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probe  data  with  a  Gaussian  distribution  whose  standard  deviation  is  related  to  the  mean 
according  to  the  rule: 

a  =  (19) 


In  Figure  8(b)  we  plot  the  relationship  between  log(/r  —  ^o)  and  log(d)  for  count  rate 
data  measured  with  the  Micro-R  probe  at  different  distances  from  the  source.  By  using  a 
linear  fit  to  this  data  we  obtain  the  following  model  for  the  mean  values  of  the  data  as  a 
function  of  the  distance  from  the  source: 


89322 


(20) 


Micro-R  count  rate  probe 


Figure  8:  The  characterisation  of  Micro-R  probe  data 


5  Probabilistic  solutions 

Because  radioactive  decay  is  an  inherently  stochastic  process,  it  is  reasonable  to  ap¬ 
ply  probabilistic  techniques  for  the  radiological  source  estimation  problem.  Apart  from 
statistical  variations  in  radiation  emanating  from  the  source,  there  is  also  a  randomly 
fluctuating  background  radiation  level  that  affects  each  measurement.  The  background 
radiation  may  comprise  contributions  due  to  cosmic  radiation,  airborne  radioactivity,  ter¬ 
restrial  radiation,  and  natural  radioactivity  of  constituent  material  of  the  detector  itself, 
other  equipment  and  shielding.  The  unknown  radiation  source  has  to  be  localised  using 
measurements  collected  in  the  presence  of  this  randomly  fluctuating  background  radiation 
level.  This  can  be  considered  as  a  parameter  estimation  problem  where  the  parameter  vec¬ 
tor  [xo  Ho  l]  has  to  be  estimated  based  on  the  set  of  measurements  {zk  ■  k  =  1, . . . ,  N}. 

In  many  parameter  estimation  problems,  the  measurements  Zk  are  modelled  as  a  non¬ 
linear  function  of  the  parameter  vector  x  with  additive  noise  as  follows  [5]: 


Zk  =  hkfx.)  +  Vk- 


(21) 
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While  this  model  may  be  appropriate  for  situations,  such  as  estimating  a  signal  con¬ 
taminated  with  noise,  it  is  somewhat  different  in  the  problem  at  hand.  In  the  radiological 
source  estimation  problem,  the  variance  in  measurements  is  largely  due  to  statistical  fluc¬ 
tuations  in  the  number  of  disintegrations  per  unit  time  in  the  source  itself,  rather  than 
due  to  the  measurement  noise.  Therefore,  when  the  above  model  is  considered  for  the 
radiological  source  estimation  problem,  /ifc(x)  may  be  considered  as  the  mean  value  of  the 
random  parameter  vector  of  interest  and  the  contribution  due  to  its  random  fluctuations 
and  noise. 

In  many  estimation  problems,  the  measurement  noise  in  Equation  (21)can  be  rea¬ 
sonably  modelled  using  a  Gaussian  density.  A  Gaussian  noise  model  is  mathematically 
attractive  and  usually  leads  to  convenient  analytical  solutions.  Popular  nonlinear  filtering 
techniques  such  as  the  extended  Kalman  filter  (EKF)  and  the  unscented  Kalman  filter 
(UKE)  are  based  on  such  modelling. 

Radiation  survey  data  typically  represent  samples  drawn  from  a  Poisson  distribution 
as  was  observed  in  the  experimental  data  collected  with  the  gamma  probe.  Since  the 
Gaussian  extends  from  — oo  to  -|-oo,  it  cannot  truly  represent  the  radiation  survey  mea¬ 
surements,  which  are  always  non-negative.  For  large  values  of  measurements,  however,  a 
Gaussian  with  a  =  is  almost  identical  to  a  Poisson  distribution  with  the  same  mean^. 
The  Gaussian  assumption  will  be  more  applicable  for  large  values  of  /ifc(x)  =  be¬ 

cause  the  truncation  point  will  be  far  from  By  approximating  the  Poisson  distributed 

^  k 

measurements  by  the  Gaussian,  we  are  able  to  model  the  data  reasonably  accurately  while 
enjoying  the  mathematical  convenience  afforded  by  the  use  of  the  Gaussian  distribution. 

Hence,  for  normalised  dose  rate  data,  we  propose  the  following  model  of  the  measure¬ 
ment  likelihood  function,  given  the  parameter  vector  x: 


p(zfc|x)  =  V [zk]  \k{y^)) 

K,  AA(2:fc;  Afe(x),  Afc(x)) 


(22) 

(23) 


where  /U,  cr^)  denotes  a  Gaussian  density  with  mean  ^  and  variance  and 


Afc(x) 


I 

{xk  -  xo)2  -h  {vk  -  yoY 


+  Afe. 


(24) 


Here  Af,  represents  the  mean  of  the  background  radiation,  that  is,  the  mean  of  the  dose-rate 
measurements  when  /  =  0  or  when  dk  =  {xk  —  xqY  +  {yk  —  yoY  is  infinite^. 

In  practice  the  measurements  are  subjected  to  a  detection  threshold.  Typically  a  de¬ 
tection  (measurement)  is  reported  only  if  it  exceeds  a  multiple  of  the  background  radiation 
mean  Xf,.  Thus  the  approximation  in  (23)  is  valid  for  all  practical  purposes.  For  this  case 
we  define  an  approximate  signal-to-noise  (SNR)  ratio  in  dB  as: 


SNRfc[(iH]  lOlogAfc  10  log  (25) 

“fc 

^The  Micro-R  probe  data  had  a  difTerent  relationship  between  a  and  /i.  The  reason  for  this  will  be 
investigated 

^Note  that  the  likelihood  function  p{zk\x)  was  modelled  by  a  truncated  Gaussian  density  in  [6]. 
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6  Cramer-Rao  analysis 

The  problem  is  to  estimate  the  parameter  vector  x  using  the  collected  measurements 
{zk  :  k  =  and  assuming  that  the  source  is  localised  within  a  specified  region  of 

interest  (this  information  will  be  our  prior  knowledge  in  the  Bayesian  sense). 


6.1  Derivation 


In  this  section  we  compute  the  theoretically  best  achievable  second  order  error  per¬ 
formance  of  any  unbiased  estimator  of  the  parameter  vector.  Index  k  is  assigned  to 
the  parameter  estimate  to  indicate  that  k  measurements  were  used  in  estimation.  This 
estimation  error  lower  bound,  known  as  the  Cramer-Rao  bound  (CRB),  will  be  a  good 
indicator  of  observability  of  location  and  activity  of  the  source. 

The  CRB,  denoted  by  C^,  provides  a  lower  bound  on  the  mean-square  error  (MSE) 
matrix  of  any  unbiased  estimator  x^.  The  MSE  matrix  is  defined  as: 

Efc  =  E  |(xfc  -  x)  (xfc  -  x)'^|  ,  (26) 

where  E  is  the  expectation  operator.  By  definition  the  CRB  is  given  by  the  inverse  of  the 
information  matrix  J^,  i.e. 

=  (27) 

Here  the  matrix  inequality  indicates  that  'Sk  —  Ck  is  a  positive  semi-definite  matrix.  In 
the  Bayesian  framework,  the  information  matrix  consists  of  twofold  contribution:  prior 
information  and  measurement  contribution. 


Assuming  Gaussian  model  in  (23),  the  information  matrix  can  be  computed  in  the 
Bayesian  framework  via  the  following  recursive  formula  [7,  8] : 

1 


Jfc  —  J k—l  T 


Afc(x; 


■Hfc(x)^Hfc(x). 


(28) 


Here  Hfc(x)  is  the  Jacobian  of  the  nonlinear  measurement  function  Afc(x)  defined  in  (24), 
that  is 


Hfc(x)  = 

9Afc(x) 

dxo 

aAfc(x)  9Afc(x) 

dyo  dl  \  ’ 

(29) 

with  individual  terms: 

5Afc(x)  _ 

2/(xo  -  Xk) 

(30) 

dxo 

[{xo 

-  XkY  +  {yo  -  ykYf 

5Afc(x)  _ 

27(2/0  -  2/fc) 

(31) 

dyo 

[(2^0 

-  xkY  +  (2/0  -  ykYf 

9Afc(x) 

1 

(32) 

dl 

(a:o  - 

xkY  +  (2/0  -  ykY' 

The  recursion  in  (28)  is  initialised  with  Jq  which  corresponds  to  our 

prior  information  (in 

the  Bayesian  sense)  about  the  parameter  vector.  Assuming  a  Gaussian  prior  distribution 
with  mean  xq  and  covariance  Pq,  we  have  that  Jq  =  Pq^-  The  prior  mean  xq  and 
covariance  Pq  may  be  obtained,  for  example,  based  on  intelligence  information.  If  prior 
information  is  weak,  Pq  will  take  large  values  which  in  the  limiting  case  yields  Jo  =  0. 
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6.2  Analysis 

Let  us  consider  first  a  scenario  where  the  measurements  are  collected  uniformly  along  a 
circular  path  with  the  source  positioned  in  the  centre  of  the  circle.  The  first  measurement 
is  taken  at  the  angle  of  45°  from  the  source  with  respect  to  the  x-axis,  with  observer 
moving  anti-clockwise.  The  activity  of  the  source  is  fixed  at  /  =  10®,  while  the  radius  of 
the  circle  takes  values  of  562.3m,  177.8m  and  56.2m,  which  correspond  to  SNR  values  of 
5dB,  15dB  and  25dB,  respectively.  The  initial  covariance  is  diagonal: 

Po  =  diag  ( [o-2  al  aj] )  (33) 

with  (Jx  =  (Jy  =  250m  and  aj  =  5  ■  10®.  The  choice  if  these  values  in  our  example  is 
arbitrary.  The  values  for  ax  and  ay  are  chosen  to  reflect  prior  knowledge  that  the  source 
is  located  with  probability  95%  within  an  elliptical  region  whose  axes  are  approximately 
2ax  and  2ay.  The  prior  information  on  source  intensity  is  weak  hence  aj  is  set  to  a  large 
value.  The  best  achievable  standard  deviation  of  the  estimation  error  in  position,  defined 
as  y^Ck[l,l]  +  Ck[2,2],  is  plotted  in  Figure  9(a)  for  the  three  values  of  SNR.  Likewise, 
Figure  9(b)  shows  the  best  achievable  standard  deviation  of  estimation  error  in  activity, 
Y^Cfc[3,3].  Observe  how  in  both  figures  the  error  standard  deviation  reduces  as  we  collect 
more  measurements.  The  importance  of  having  a  high  enough  SNR  is  also  evident:  after 
processing  all  60  measurements,  at  SNR=15  (versus  5dB)  we  can  achieve  the  standard 
deviation  of  positional  estimation  error  of  5m  (versus  50m).  CRB  analysis  allows  us  to 
quantify  the  best  achievable  estimation  accuracy  for  different  values  of  SNR  and  for  a 
given  scenario. 


7  Estimation  algorithms  and  their  performance 

7.1  Algorithms 

In  this  section  we  describe  several  estimation  algorithms  for  estimation  of  parameter 
vector  X.  Parts  of  this  section  have  been  reported  in  [9]. 

7.1.1  Maximum  likelihood  estimation 

Since  we  deal  here  with  parameter  estimation,  the  first  algorithm  to  consider  is  the 
maximum  likelihood  estimator  (MLE).  The  maximum  likelihood  estimator  (MLE)  is  a 
block  algorithm  as  it  operates  on  the  accumulated  set  of  all  previous  measurements.  It 
does  not  use  any  prior  information  (hence  referred  as  non-Bayesian)  but  is  widely  used 
for  parameter  estimation  because  if  an  asymptotically  unbiased  and  minimum  variance 
estimator  exists  for  large  sample  sizes,  it  is  guaranteed  to  be  the  MLE  [10].  The  MLE  is 
determined  as  the  vector  which  maximises  the  likelihood  function  p{zi, . . . ,  2;Ar|x). 

=  argmax  p{zi, . . . ,  ZArjx).  (34) 

X 

Assuming  that  measurements  are  mutually  independent,  with  the  likelihood  of  each 
measurement  given  by  Equation  (23),  the  likelihood  function  becomes: 
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Number  of  measurements 


(a) 


Figure  9:  CRB  analysis:  the  error  standard  deviation  in  (a)  position  and  (h)  aetivity 


N 


p{zi, . . .  ,ZAr|x)  = 


\/27rAfc(x) 


:  exp 


1  {zk  -  M(x)y 

2  Afc(x) 


By  taking  the  logarithm  of  both  sides  of  Equation  (35)  we  obtain: 


logp(zi,  ZAr|x)  =  (  -^  )  X]  (27rAfc(x))  + 


N 


k=l 


(zk  -  \k{^)y 
Afc(x) 


Then  the  MLE  of  Equation  (34)  becomes: 

{zk  -  Xk{'x-)Y 


=  arg  min  |  ^ 

^  k=l 


Afc(x) 


+ 


iV 

^log(27rAA:(x))|. 


k=l 


(35) 


(36) 


Minimisation  in  (34)  can  be  performed  by  numerical  methods  [11];  our  implementation 
is  based  on  MATLAB®  built-in  routine  fminsearch.  Being  a  non-Bayesian  estimator,  the 
MLE  is  expected  to  reach  the  Bayesian  CRB  (discussed  in  Section  6)  only  for  large  enough 
sample  size  (large  A^)  when  the  influence  of  prior  knowledge  has  diminished. 
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7.1.2  Least  squares  (LS)  with  reduced  search  dimension 


The  least  squares  (LS)  method  is  another  commonly  used  estimation  procedure.  If  Zk 
is  the  A:th  measurement  and  /ifc(x)  its  prediction,  the  LS  estimate  is  defined  as: 

N 

x^"^  =  argmin  Y]  (zfc  - /rfc(x))^  .  (37) 

k=l 


Equation  (37)  does  not  make  any  assumptions  about  the  measurement  errors  [5].  If 
these  errors  are  independent  and  identically  distributed  zero-mean  Gaussian  random  vari¬ 
ables,  then  the  LS  estimate  coincides  with  the  ML  estimate. 


Typically,  one  may  directly  minimise  the  sum  of  squared  errors  by  searching  in  the 
three  dimensional  space  spanned  by  (xo,yO)7))  for  example  using  the  fminsearch  routine. 
We  consider  an  alternative  approach  developed  by  Takumi  et  al.[12]  in  which  the  search 
is  carried  out  in  only  the  two  dimensional  (xo,yo)  space. 

Let  /  denote  the  sum  of  the  squared  error  terms. 

TV  /  j^2 

f{xo,  yo,I)  =  '^[zk-  (38) 

k=i  ^ 


where  (x^  —  xq)^  -|-  {uk  —  VoY  is  the  radial  distance 

to  the  putative  radiation  source.  At  the  optimal  value  of 
estimation  error  /, 


from  the 


measurement  point 
that  minimises  the 

(39) 


From  this  equation,  the  optimal  value  of  I  can  be  obtained 


as  a  function  of  (xo,yo))  i-e. 


I  = 


2^k=l 


Zk 

if 


1 

2-^k=l  ^ 

k 


(40) 


Equation  (40)  can  be  used  to  compute  the  value  of  I  at  all  possible  (xo,yo)  coordinates 
within  an  area  where  the  unknown  source  is  assumed  to  be  located.  By  substituting  the 
computed  I  in  (38),  the  estimation  error  is  computed  at  each  point.  The  (xo,2/o)  value 
pair  corresponding  to  the  minimum  value  of  /  is  chosen  as  the  most  likely  source  location 
within  the  search  space.  Because  this  approach  reduces  the  search  space  dimension  from 
3D  to  2D,  the  minimum  can  be  obtained  relatively  rapidly.  The  search  may  be  carried 
out  in  more  than  one  stage,  where  a  large  area  is  first  searched  at  a  coarse  resolution  to 
find  the  most  likely  subarea,  which  will  then  be  searched  using  a  finer  resolution  search. 
Instead  of  simply  accepting  the  minimum  solution  returned  by  the  algorithm,  it  may  be 
informative  to  look  at  a  plot  of  the  squared  error  surface,  which  could  show  other  likely 
source  positions.  A  good  example  is  when  data  is  collected  along  a  linear  trajectory.  While 
the  algorithm  may  return  the  location  corresponding  to  the  minimum  error,  a  plot  of  the 
error  surface  may  show  that  an  image  solution  also  exists. 


Since  this  is  a  nonlinear  LS  estimation  problem  where  according  to  Equations  (23)  and 
(24),  additive  noise  has  a  non-constant  variance,  a  theoretical  prediction  of  the  estimation 
accuracy  would  be  very  difficult  to  derive. 
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7.1.3  An  approximate  recursive  LS  algorithm 

Although  the  least  squares  approach  is  capable  of  localising  the  source  accurately,  being 
a  batch  algorithm,  it  is  unattractive  for  real-time  operational  use.  A  batch  algorithm  re¬ 
quires  increasing  amounts  of  memory  and  processing  time  as  the  number  of  measurements 
is  increased.  A  recursive  algorithm,  in  contrast,  requires  essentially  the  same  amount  of 
processing  time  and  memory  for  each  new  measurement  irrespective  of  the  total  number 
of  measurements  and  is  attractive  for  real-time  applications.  While  a  recursive  solution 
is  straightforward  and  widely  used  in  linear  least  squares  problems,  it  is  not  the  case  for 
non-linear  least  squares  problems.  Some  possible  approaches  to  derive  recursive  non-linear 
least  squares  solutions  are  described  in  the  literature  [13,  14,  15]. 

While  we  intend  to  investigate  these  approaches  in  future  and  apply  them  to  the 
radiation  source  estimation  problem,  here  we  describe  a  simple  inexact  recursive  least 
squares  solution  that  we  developed  for  initial  experimentation.  Although  it  is  only  an 
approximate  solution,  it  was  able  to  provide  reasonably  accurate  source  estimates  when 
applied  to  the  set  of  real  measurements  currently  available  to  us.  These  results  will  be 
presented  in  Subsection  8.5.  To  derive  this  approximate  recursive  solution,  we  rewrite 
Equation  (40)  as  follows: 


In 


l^k=l 


Zk 

'7^ 


1 

l^k=l  ^ 

k 


A/jv 

Vn 


(41) 


where  In  is  the  I  estimate  after  observing  the  Nth.  observation,  zn,  and  Nn  and  Vn  are 
the  numerator  and  the  denominator  of  Equation  (41).  Therefore,  the  estimate  of  I  after 
the  {N  +  l)th  observation  can  be  expressed  as: 


where: 


^AT+l 


A/jv+1 

Vn+1 


Mn+1  =  Mn  + 


^jV+l 

p2 

-^TV+l 


Vn+1  =Vn  +  ^4 — 

^N+l 


(42) 


(43) 

(44) 


Unfortunately,  Equation  (38)  cannot  be  transformed  into  a  recursive  form.  Therefore, 
we  tried  the  following  recursive  form,  which,  while  not  exactly  equivalent  to  Equation  (38), 
tries  to  approximate  it  to  a  limited  extent: 


Jn+i 


In  +  ~ 


Av±]_  \ 

Rl+i 


(45) 


The  source  estimates  after  observing  the  (A-|-l)th  measurement  are  computed  by  searching 
for  the  minimum  of  fN+i-  Rigourous  approaches  to  recursive  nonlinear  least  squares  will 
be  investigated  in  future  work. 
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7.1.4  EKF  and  UKF 

The  radiological  source  localisation  problem  can  be  formulated  within  the  Bayesian 
filtering  framework  [8] .  Since  the  parameter  vector  is  fixed,  the  state  dynamic  equation  is 
trivial.  The  measurement  equation  (Eqn  23)  is  highly  nonlinear  and,  therefore,  we  need 
to  consider  a  nonlinear  filtering  solution. 

The  goal  of  filtering  (the  sequential  Bayesian  estimation)  is  to  determine  the  posterior 
pdf  p{x\zi, . . . ,  Zfc),  for  k  =  1,2,....  Both  the  extended  Kalman  filter  (EKF)  [5]  and  the 
unscented  Kalman  filter  (UKF)  [16]  approximate  this  posterior  pdf  by  a  Gaussian  density 
The  mean  and  covariance  of  this  Gaussian  pdf  are  computed  recursively: 
in  the  case  of  the  EKF  via  the  local  linearisation  of  the  measurement  function  /ifc(x);  in 
the  case  of  the  UKF  via  an  approximation  of  the  posterior  pdf  by  a  set  of  deterministically 
chosen  sample  (or  sigma)  points. 

Both  EKF  and  UKF  are  based  on  approximations  and  hence  suboptimal  and  prone  to 
divergences  (as  we  will  see  in  Sec. 7. 2).  However,  the  UKF  is  a  more  accurate  algorithm 
since  it  captures  the  mean  and  covariance  of  the  underlying  posterior  pdf  exactly  up  to 
the  second  order  of  nonlinearity.  Both  EKF  and  UKF  implementations  contain  a  small 
and  equal  amount  of  process  noise. 

7.1.5  EKF  equations 


Xfc|fc-i  —  ffc-i(xfc_i|fc_i)  (46) 

Pfc|fc-i  =  Qfc-i  +  Ffc_iPfc_i|fc_iF^_;^  (47) 

are  the  state  prediction  and  state  prediction  covariance  based  on  measurements  up  to  the 
{k  —  l)th  measurement.  The  notation  k\k  —  1  denotes  the  prediction  of  the  kth  estimate 
after  observing  up  to  the  {k  —  l)th  measurement.  The  updated  estimate  after  observing 
the  kth  measurement  is  denoted  by  k\k.  Qfc_i  is  the  covariance  of  process  noise. 


^k\k  —  ^k\k-l  +  ^ki'^k  —  ^k{^k\k-l)) 

(48) 

^ k\k  Ffc|fc— 1  “ 

(49) 

are  the  updated  state  and  state  covariance  after  receiving  the  kth  measurement  z^- 

Sfc  =  HfcPfc|fc_iHfc  +  Rfc 

(50) 

is  the  innovation 

covariance,  where  the  innovation  is  defined  as: 

=  Zk  -  Ak(xk|k-i) 

(51) 

and 

Kfc  =  Pfc|fc_iH^S^^ 

(52) 

is  the  filter  gain  or  the  Kalman  gain.  is  the  covariance  of  measurement  noise. 

.1  and 

Hfc  are  the  local  liberalisations  of  nonlinear  functions  ffc_i(x)  and  Afc(x),  respectively  [8]. 


18 


DSTO-TR~1988 


Because  we  are  trying  to  estimate  the  position  and  the  strength  of  a  static  source,  the 
state  transition  equation  is  simply: 


and 


Xfc  =  Xfc_i. 


Ffc-i 


Xfc-1— 


■ 

r  a  1 

axo 

a^o 

.  ai  _ 

T 


[a;o  yo  I] 


Xfc-i— 


I. 


(53) 

(54) 

(55) 

(56) 

(57) 


Using  the  measurement  equation  (Eqn  (24))  in 

Hfc  =  [Vxj,  (Xfc)]  |xfc=x;j,n,_i 


we  get: 


H 


k  = 


a  ' 

5! 

t’ 

dl 


{x-xo)^+{y-yo)^ 


— ^k\k  —  l 


2I(x-xo)  2I{y-yo)  _  1 _ 

[(x-xo)2+(j/-j/o)^]^  [{^-xo)‘^+{y-yo)'^]^  (*-^o)^+(j/-j/o)^ 


— 1  ’ 


(58) 


(59) 

(60) 


7.1.6  UKF  equations 

In  the  case  of  the  UKF,  the  posterior  density  at  time  A:  —  1  is  assumed  to  be  Gaussian 
and  this  density  is  represented  by  a  set  of  M  sample  points  Xk-i  corresponding  weights 
Wl_^,  which  can  be  computed  according  to  the  unscented  transform  [8]. 

Then  the  predicted  state  and  the  state  prediction  covariance  are  given  by: 


M-l 

Xfc|fc-1  =  Wl_i  ■  h-i{xl-i) 

i=0 


(61) 


M-l 

Pfc|fc-1  =  Qfc-1  +  Y  [ffc-i(Xfc-i)  -Xfc|fc-i]  [ffc-i(Xfc-i)  -^k\k-i]'^  ■  (62) 

i=0 

The  predicted  measurement  is  computed  as: 

M-l 

^k\k-i  =  Y  n-Mxi\k-i)-  (63) 

i=0 
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Upon  receiving  the  kth  measurement  z^,  the  state  and  the  state  covariance  are  updated 
as  follows: 


^k\k  —  ^k\k-l  +  ^ki^k  —  ^k\k-l) 

(64) 

P  k\k  1 

(65) 

where 

Kfc  = 

(66) 

is  the  filter  gain  and 

Sfc  —  rtfc  ~kPzz 

(67) 

is  the  innovation  covariance.  Y^xz  and  'Pzz  are  computed  using: 

M-l 

Pxz  = 

tF^-i(Xfc|fc-i  -ifc|fc-i)(Afc(xl|fc-i)  -  h\k-if 

i=0 

M-l 

(68) 

Pzz  = 

Y  ^Li(hfc(Xfc|fc_i)  -  Zfc|fc-i)(Afc(xl|fc_i)  -  h\k-if- 

i=0 

(69) 

7.2  Simulation  results 

In  order  to  compare  the  performance  of  proposed  algorithms,  we  carried  out  100  Monte 
Carlo  runs  using  the  scenario  shown  in  Figure  10,  where  the  crosses  indicate  the  locations 
where  the  measurements  are  taken,  and  the  star  marks  the  location  of  the  source.  The 
source  activity  is  set  to  /  =  10®,  so  that  at  the  closest  point  of  approach  (205m  away  from 
the  source),  the  SNR  is  14dB.  Other  parameters  used  in  simulations  are:  cr^  =  (Jy  =  250m, 
aj  =  10®,  N  =  60,  Afe  =  0.85.  The  MLE,  being  a  block  algorithm,  needs  to  be  re-run  every 
time  a  new  measurement  becomes  available,  in  order  to  show  the  estimation  error  versus 
time.  In  addition,  it  is  required  to  accumulate  at  least  33  measurements  to  run  the  MLE 
(otherwise  it  becomes  unstable). 

Eigures  11  and  12  show  the  standard  deviation  of  the  estimation  error  for  source 
position  and  activity,  respectively.  Analysis  of  the  error  performance  results  leads  to  the 
following  observations.  The  EKE  is  diverging  and  is  considered  as  unsuitable  for  this 
application.  The  UKE  follows  the  trend  of  the  CRB,  but  is  somewhat  worse  than  the 
bound.  The  MLE  is  indeed  an  asymptotically  efficient  (unbiased  and  minimum  variance) 
estimator:  its  bias  is  approaching  zero  (not  shown  here)  and  its  error  standard  deviation 
is  approaching  the  CRB  as  more  measurements  become  available  for  estimation.  Finally 
we  point  out  that  the  theoretical  CRB  is  indeed  a  good  predictor  of  the  best  achievable 
error  performance. 


8  Application  to  real  radiological  survey  data 

8.1  Holsworthy  field  trial  data 


20 


To  verify  the  theoretical  analyses  we  applied  them  to  a  limited  amount  of  real  radi¬ 
ological  survey  data  collected  at  Holsworthy  Barracks,  NSW  in  May  2005  during  a  field 
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Figure  1 0:  The  scenario  used  for  simulation 
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Figure  1 1 :  The  estimation  error  in  position 
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Figure  12:  The  estimation  error  in  source  activity 
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trial  conducted  to  evaluate  the  LCAARS  system. 

Two  point  gamma  radiation  sources;  a  ®®Co  source  (source  reference  number  14) 
of  approximately  0.3GBq  activity  emitting  1.17  and  1.33  MeV  gamma  radiation  and  a 
^^^Cs  source  (source  reference  number  23)  of  approximately  26.8  GBq  emitting  0.662  MeV 
gamma  radiation  were  used  in  the  field  trial.  One  set  of  dose  rate  data  and  one  set  of 
count  rate  data  were  collected  with  each  one  of  the  two  radioactive  sources  emplaced.  One 
set  of  background  data  was  also  collected  with  each  detector  in  the  absence  of  the  sources. 
In  addition  to  the  dose  rate  or  count  rate  data,  GPS  receiver  readings  of  longitude  and 
latitude  of  the  detector  position  were  also  recorded.  Because  all  estimation  algorithms  as¬ 
sume  local  Gartesian  coordinates,  longitude  and  latitude  data  were  converted  to  Gartesian 
coordinates  using  the  M_MAP  mapping  package.  Unfortunately,  due  to  some  technical 
problems,  the  GPS  coordinates  of  the  two  sources  had  not  been  recorded  accurately  and 
are  known  only  approximately. 


8.2  Deterministic  solutions 


The  inverse  distance  square  law-based  analytical  solution  developed  in  Section  3  was 
applied  to  the  real  data  from  the  Holsworthy  trial.  The  four-point  algorithm  described 
therein  requires  as  input  the  x  and  y  coordinates  at  four  arbitrary  points  and  the  radiation 
survey  data  measured  at  these  four  positions.  For  the  purpose  of  illustration,  the  x  and  y 
coordinates  and  data  at  four  points  selected  from  the  ^^^Gs  count  rate  data  set  are  shown 
in  Table  3. 


Table  3:  An  example  set  of  four  measurements  from  the  Cs  eount  rate  data  used  to 
demonstrate  souree  estimation  using  the  four  point  estimation  algorithm 


xi=1946.1 

2/1=396.1 

T>i=3950 

X2=2031.0 

2/2=627.0 

112=2820 

X3=1949.6 

2/3=478.6 

113=8860 

X4=1854.9 

2/4=426.8 

114=2150 

Based  on  these  four  data  points  the  algorithm  estimated  the  source  position  as  (2004.3, 
491.2),  which  is  in  very  good  agreement  with  the  true  source  position.  In  Figure  13  the 
blue  circles  show  all  positions  at  which  measurements  were  collected,  with  the  diameter  of 
each  circle  indicating  the  relative  magnitude  of  the  corresponding  measurement.  The  four 
red  squares  show  the  four  measurements  described  in  Table  3.  The  green  dot  and  the  black 
triangle  denote  the  true  source  location  and  the  estimated  source  location,  respectively. 

To  test  the  reliability  of  the  algorithm,  a  collection  of  30  four-point  groups  were  chosen 
from  the  ^^'^Gs  count  rate  data  and  each  group  of  four  points  was  successively  input  to 
the  algorithm.  The  estimates  of  source  location  obtained  in  these  trials  are  shown  in 
Figure  14.  Except  for  one  outlier  (  at  (4205,  229))  out  of  the  thirty  estimates,  the  others 
are  clustered  closely  around  the  true  source  position.  It  is  remarkable  that  each  estimate 
is  based  on  only  four  measurements.  Experimentation  with  the  data  shows  that  estimates 
may  become  unreliable  if  two  or  more  of  the  four  measurement  points  are  too  close  to  each 
other  which  is  the  case  in  the  outlier  in  Eigure  14. 
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Figure  13:  An  example  application  of  the  four  point  based  deterministic  source  estimation 
algorithm. 
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The  estimation  algorithm  used  here  was  developed  assuming  the  measurements  to 
be  deterministic.  Because  the  radioactive  decay  is  a  stochastic  process  rather  than  de¬ 
terministic,  the  inverse  square  law  is,  actually,  only  applicable  to  the  average  values  of 
measurements,  rather  than  to  the  individual  measurements.  Therefore,  we  can  expect  the 
algorithm  to  produce  more  accurate  and  less  variable  source  estimates  if  the  mean  values 
of  measurements  at  the  sampling  points  are  used.  We  cannot  investigate  this  hypothesis 
with  the  raw  Holsworthy  trial  data,  as  only  a  single  measurement  had  been  collected  at 
most  sensor  locations.  During  the  field  trial  planned  for  October  2006,  we  intend  to  collect 
multiple  measurements  at  each  measurement  location  which  would  allow  us  to  obtain  the 
mean  value  at  each  position. 


8.3  EKF,  UKF  and  MLF 


In  subsection  4.1  it  was  shown  that  the  normalised  dose  rate  data  were  well  modelled 
by  a  Poisson  distribution.  An  appropriate  model  for  measurement  likelihood  was  proposed 
in  Equation  24. 

In  order  to  apply  the  postulated  model  to  real  data  collected  in  the  field  trial  we 
require  the  mean  value  \h  corresponding  to  the  background  data  measured  in  situ.  The 
normalised  background  dose  rate  data  collected  with  the  gamma  probe  in  the  field  trial 
and  a  Poisson  pdf  fitted  to  this  data  are  shown  in  Figure  15.  From  these  data  Xh  was 
estimated  to  be  0.8376,  which  is  slightly  larger  than  the  value  of  0.7946  estimated  from 
laboratory  data.  The  shielding  provided  by  the  concrete  walls  and  the  top  floors  of  the 
building  may  be  the  reason  for  the  lower  background  level  inside  the  laboratory. 


Holsworthy  Trial  Gamma  probe  background  data 


Figure  15:  The  normalised  histogram  for  baekground  dose  rate  eolleeted  at  Holsworthy 
army  barraeks  and  the  Poisson  distribution  fitting 

Before  we  applied  the  estimation  algorithms  (MLF  and  UKF)  we  have  discarded  all 
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dose  rate  measurements  with  a  normalised  value  lower  than  three.  Both  estimation  al¬ 
gorithms  were  given  the  same  prior  information:  the  region  where  the  source  is  located 
(500m  X  500m)  and  the  guessed  value  of  source  activity.  For  the  MLE  this  prior  in¬ 
formation  was  used  to  initialise  the  fminsearch  routine,  while  for  the  UKF  this  was  the 
initial  value  of  the  state  vector.  In  addition,  the  elements  of  Pq  (required  for  UKF)  were: 
(^x  =  CTy  =  125m,  aj  =  10^.  The  results  of  our  analysis  are  shown  in  Figure  16. 


X[m] 

Figure  16:  Trial  data  analysis:  Measurement  loeations/intensities  (o),  true  (*)  and 
estimated  (D  and  V )  souree  locations 

Blue  circles  in  Figure  16  show  the  locations  and  intensities  of  measured  dose- rates. 
The  centre  of  each  circle  indicated  the  location  {xk,yk),  while  the  radius  of  each  circle 
corresponds  to  z^-  The  approximate  true  location  of  the  source  is  marked  with  a  green  * 
symbol,  the  MLE  and  UKE  estimates  with  a  red  triangle  and  a  black  square,  respectively. 
Because  the  EKE  diverged  and  failed  to  return  a  meaningful  estimate  these  results  are 
not  shown  here. 

Eigure  17  shows  the  actual  measurements  Zk  (circles)  versus  Afc(x)  (crosses),  where  x 
is  the  MLE  estimate.  The  latter  represents  the  mean  of  p(zfc|x),  (Eqn  23).  We  observe  a 
good  agreement  between  the  measurements  and  their  prediction  based  on  the  MLE. 


8.4  LS  with  reduced  search  dimension 

We  applied  the  LS  algorithm  with  2D  search  to  the  available  radiation  survey  data. 
Eigure  18  shows  the  sum  squared  error  surface  obtained  for  ^^^Cs  count  rate  data.  The 
white  circles  show  the  locations  where  data  were  collected.  The  radius  of  each  circle  is  used 
as  a  rough  indicator  of  the  relative  magnitude  of  the  measurements.  Dark  blue  regions 
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Figure  17:  Measured  dose  rates  and  Afc(x)  eorresponding  to  the  MLE 

indicate  source  locations  that  lead  to  small  values  of  sum  squared  error  and,  therefore, 
highly  likely  to  be  the  true  source  position.  Dark  red  regions  correspond  to  positions 
that  would  lead  to  large  values  of  sum  squared  error  if  the  source  were  located  there 
and,  therefore,  highly  unlikely  to  be  the  true  source  position.  The  white  square  shows  the 
location  of  the  minimum  error,  which  is  chosen  as  the  estimated  source  location.  The  green 
asterisk  is  the  location  where  the  true  source  was  believed  to  be  located.  The  estimated 
source  position  (2000,  510),  is  in  good  agreement  with  the  true  location,  (2000,  500). 

The  distance  from  the  estimated  source  position  to  each  measurement  point  and  the  es¬ 
timated  source  activity  were  used  to  predict  the  measurement  values  at  these  measurement 
locations.  Figure  19  compares  the  measured  data  with  the  predictions. 

The  value  of  I  was  estimated  to  be  4.71  x  10^.  Using  Equation  (6)  together  with  the 
energy  response  value  of  12200  for  the  Micro-R  probe  for  ^^^Cs  radiation  the  activity  of 
the  source  is  estimated  to  be: 

M  =  (4.71  X  10^  X  6)  /  (12200  x  0.662) 

=  35  GBq 

which  is  of  the  same  order  as  the  true  source  activity  of  26.8  GBq. 

When  the  LS  algorithm  was  applied  to  dose  rate  data  collected  with  the  ^^^Cs  source 
in  place,  it  resulted  in  the  error  surface  shown  in  Figure  20.  Although  the  estimated  source 
position  (1957,  508)  is  about  43  m  away  from  the  true  source  position,  the  error  surface 
correctly  recognised  the  true  source  position  also  as  highly  likely,  as  indicated  by  the  dark 
blue  shade  corresponding  to  very  low  error  in  this  area.  How  the  measurements  predicted 
using  the  estimated  source  compare  with  the  actual  measurements  is  shown  in  Figure  21. 

Figure  22  shows  the  error  surface  obtained  for  count  rate  data  measured  with  the 
®^Co  source  in  place.  In  this  case  the  source  location  estimated  by  the  LS  algorithm 
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Figure  18:  The  squared  error  surface  obtained  for  Cs  count  rate  data.  Measurement 
points  (o)  and  the  estimated  (O)  and  true  (^)  source  positions  are  also  shown. 


Count  rate  measurements  vs  data  predicted  with  estimated  Cs137  source 


Figure  19:  Comparison  of  Cs  count  rate  measurements  and  prediction  based  on  LS 
estimation 
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Figure  20:  The  squared  error  surfaee  obtained  for  dose  rate  data.  Measurement 

points  (o)  and  the  estimated  (Ti)  and  true  (*)  souree  positions  are  also  shown. 


Figure  21: 

estimation 


,x  10 


i-3  Measurement  vs  prediction:  CS137  dose  rate  probe 

■  Measured 

■  Pred  w/estimated  source 


0  50  100  150  200  250  300  350  400 

Measurement  number 


Comparison  of  dose  rate  measurements  and  predietion  based  on  LS 


28 


DSTO-TR-1988 


is  different  to  the  true  source  position.  However,  the  dark  blue  region  near  the  true 
source  is  an  indication  that  this  area  is  also  highly  likely  to  have  the  source.  Because 
the  measurements  in  this  case  are  almost  along  a  straight  line,  an  image  solution  is,  in 
fact,  possible  due  to  non-uniqueness.  Therefore,  if  this  algorithm  is  used  in  an  operational 
situation,  it  may  be  prudent  to  search  all  areas  that  show  low  sum  squared  errors  (dark 
blue  areas  in  the  plots).  The  measured  and  predicted  data  are  compared  in  Figure  23. 
The  error  surface  corresponding  to  dose  rate  data  measured  with  the  ®°Co  source  in  place 
is  shown  in  Figure  24.  In  this  case  the  estimate  is  closer  to  the  supposed  true  source 
position. 
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Figure  22:  The  squared  error  surfaee  obtained  for  Co  eount  rate  data.  Measurement 
points  (o)  and  the  estimated  (DJ  and  true  (*)  souree  positions  are  also  shown. 

Figure  26  shows  the  error  surface  obtained  for  background  count  rate  data  measured 
in  the  absence  of  sources.  The  white  square  that  denotes  the  estimated  source  location 
is  at  the  top  right  corner  of  the  picture,  which  is  over  2  km  away  from  the  measurement 
location.  The  algorithm  correctly  recognised  in  this  case  that  it  would  not  be  possible  for 
a  source  to  be  present  close  to  the  measurement  points. 

These  results  indicate  the  LS  algorithm  to  be  a  useful  tool  to  locate  an  unknown  single 
point  gamma  radiation  source.  It  is  particularly  informative  to  look  at  the  complete  error 
surface  and  visualise  all  highly  likely  source  locations  than  to  simply  rely  on  the  single 
minimum  returned  by  the  algorithm.  This  fact  was  clearly  demonstrated  in  some  of  the 
results  presented  above.  While  the  results  of  applying  the  LS  algorithm  to  Holsworthy 
trial  data  generally  demonstrate  its  utility  in  locating  point  radiation  sources,  it  is  not 
possible  to  provide  an  exact  measure  of  the  accuracy  of  the  estimates  because  the  true 
source  locations  are  known  only  approximately. 
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Count  rate  measurements  vs  data  predicted  with  estimated  Co60  source 


Figure  23:  Comparison  of^^Co  count  rate  measurements  and  prediction  based  on  LS 
estimation 


f=l(2.-Dj|/Rp  CO60  dose  rate  probe 


2760 


Figure  24:  The  squared  error  surface  obtained  for  ^^Co  dose  rate  data.  Measurement 
points  (o)  and  the  estimated  (tHj  and  true  (*)  source  positions  are  also  shown. 
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^g-3  Measurements  vs  predicted  data:  CO60  dose  rate  probe 


Figure  25:  Comparison  of  ^^Co  dose  rate  measurements  and  predietion  based  on  LS 
estimation 


8.5  The  approximate  recursive  LS  algorithm 

An  inexact  LS  algorithm  to  estimate  the  unknown  radiation  source  based  on  radiation 
survey  data  in  a  recursive  manner  was  described  in  Subsection  7.1.3.  To  check  whether  this 
inexact  recursive  algorithm  could  produce  reasonable  source  estimates,  it  was  applied  to 
real  data.  Measured  count  rate/dose  rate  data  were  used  to  recursively  update  estimates 
of  xq,  yo  and  I  and  the  position  estimation  error  was  computed  after  each  update.  To 
compare  the  source  estimates  of  the  approximate  recursive  LS  algorithm  with  those  of  the 
exact  batch  algorithm,  the  latter  algorithm  was  also  run  repeatedly,  varying  the  size  of 
the  batch  from  one  to  the  full  set  of  measurements.  Plots  comparing  the  progression  of 
estimates  for  the  approximate  recursive  algorithm  and  the  exact  batch  algorithm  for  the 
case  of  ®^Co  count  rate  data  are  shown  in  Figures  27(a)-(c).  The  corresponding  progression 
of  position  errors  for  the  two  algorithms  are  compared  in  Figure  27(d).  These  plots  indicate 
that  the  estimates  of  the  recursive  algorithm  match  those  of  the  batch  algorithm  for  these 
data.  Similar  agreement  was  obtained  for  the  other  data  sets. 

The  estimates  fluctuate  rapidly  until  about  the  70th  measurement  in  the  data  shown, 
after  which  the  estimates  stabilise  to  the  final  value.  From  Figure  19  we  see  that  this 
corresponds  to  a  strong  measurement  peak.  The  early  measurements  are  mostly  low  count 
rate  data  within  the  background  level.  Furthermore,  because  MLE  is  only  asymptotically 
unbiased,  estimates  obtained  with  a  small  number  of  data  points  may  be  biased.  The 
recursive  algorithm  seems  to  produce  more  fluctuations  than  the  batch  algorithm,  initially, 
but  after  processing  the  first  strong  measurement  peak  it  also  converges  rapidly  to  the 
correct  estimates. 


The  recursive  algorithm  is  quite  attractive,  as  it  can  be  used  in  real  time.  In  a  real 
operational  scenario,  obviously,  the  position  error  is  not  known.  However,  the  plots  of 
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Figure  26:  The  squared  error  surface  obtained  for  background  measurements.  Measure¬ 
ment  points  are  denoted  by  (o)  and  the  estimated  source  position  by  (fT) 

estimates  xq,  yo  and  /,  similar  to  those  in  Figures  27(a)-(c)  can  be  generated  in  real  time. 
If  the  estimates  stop  fluctuating  and  stabilise  to  some  values,  this  may  be  considered  as 
an  indication  of  convergence  to  a  reliable  source  estimate.  Once  this  is  observed,  it  may 
be  useful  to  move  around  this  estimated  source  position  and  collect  further  measurements 
to  refine  the  estimate.  If  the  estimated  source  position  were  correct,  moving  towards  this 
estimated  location  should  produce  measurements  that  are  stronger  on  average.  Of  course, 
moving  towards  the  estimated  source  should  be  done  only  as  long  as  the  measured  count 
or  dose  rates  do  not  exceed  any  predetermined  maximum  safe  exposure  levels. 

Although  this  recursive  algorithm  performed  well  with  the  limited  real  data  available 
from  the  Holsworthy  trial,  because  it  is  an  inexact  algorithm,  much  more  testing  with  real 
data  should  be  carried  out  to  validate  it. 


9  Conclusions 

Several  algorithms  for  estimating  a  single  radiological  point  source  were  developed  and 
studied  using  simulated  and  real  radiological  survey  data. 

A  deterministic  analytical  solution  based  on  the  assumption  of  the  inverse  square  law 
was  able  to  provide  rough  source  estimates  using  measurements  collected  at  just  four 
arbitrary  points.  Because  the  inverse  square  law  is  applicable  to  a  radiological  source  only 
in  an  average  sense,  this  algorithm  can  be  expected  to  produce  more  stable  and  accurate 
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Progression  of  estimate  (Cs137  all)  Progression  of  estimate  (Cs137  all) 


Measurement  number  Measurement  number 


(a)  The  xq  estimate 


(b)  The  yo  estimate 


Progression  of  estimate  (Csl  37  all)  g^ror  tor  Csl  37  count  rate  data 


(c)  The  I  estimate 


(d)  The  position  error 


Figure  27:  The  comparison  of  the  batch  LS  algorithm  with  the  approximate  recursive  LS 
algorithm  for  the  case  of  Cs  count  rate  data 
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estimates  if  averages  of  multiple  sensor  measurements  collected  at  each  position  rather 
than  individual  readings  are  used  as  input  to  the  algorithm. 

The  theoretical  estimation  lower  bound  known  as  the  Cramer-Rao  bound  (CRB)  was 
computed  to  visualise  the  observability  of  the  source  location  and  the  activity.  The  max¬ 
imum  likelihood  estimate  (MLE)  was  obtained  and,  using  simulation,  it  was  shown  to 
approach  the  theoretical  bound  predicted  by  the  CRB  analysis.  When  applied  to  real 
measurement  data,  the  MLE  was  in  good  agreement  with  the  known  truth  values.  Al¬ 
though  the  MLE  is  an  asymptotically  optimal  algorithm  that  approaches  the  theoretical 
bound  when  the  number  of  measurements  is  large,  being  a  batch  algorithm,  it  is  less 
attractive  for  operational  use. 

The  unscented  Kalman  filter  (UKE)  and  the  extended  Kalman  filter  (EKE)  were  also 
studied.  While  the  estimation  error  of  the  UKE  was  worse  than  that  of  the  MLE,  it 
produced  estimates  that  were  reasonably  close  to  the  true  source  values.  The  EKE  diverged 
and  failed  to  produce  acceptable  estimates. 

A  least  squares  approach  that  uses  an  exhaustive  search  in  a  two  dimension  search  space 
was  investigated  and  shown  to  return  good  estimates.  An  inexact  recursive  version  of  this 
approach  was  developed.  While  this  recursive  algorithm  produced  reasonably  accurate 
source  estimates  when  applied  to  the  limited  set  of  real  data  available,  its  validity  needs 
to  be  verified  with  more  real  data. 

The  radiological  source  localisation  problem  is  a  complex  problem.  The  goal  of  re¬ 
search  in  this  area  is  to  be  able  to  reliably  estimate  multiple  possibly  moving  sources  and, 
ultimately,  track  a  cloud  of  radioactive  particles.  This  report  described  our  initial  work 
towards  understanding  this  difficult  problem  by  first  studying  the  relatively  simple  case  of 
a  single  fixed  radiological  source. 

We  list  below  some  possible  research  directions  for  future  work. 

1.  Joint  detection  and  localisation  of  multiple  moving  point  sources 

Current  work  was  based  on  prior  knowledge  that  a  single  static  point  source  is 
active  in  a  region  of  interest.  In  reality  this  may  not  be  realistic,  and  we  may 
need  to  estimate  both:  (a)  the  number  of  point  sources  and  (b)  their  location  (with 
possibly  other  attributes).  The  sources  may  also  be  moving.  Einally,  it  might  be 
of  interest  to  investigate  the  influence  of  the  measurement  location  uncertainty  (e.g. 
GPS  accuracy  is  about  3m)  on  source  localisation. 

2.  Localisation  of  point  sources  in  a  non-homogeneous  radiation  environment 

Current  work  assumes  that  the  attenuation  of  the  radiation  is  equal  in  all  directions 
(e.g.  open  space  homogeneous  environment).  In  urban  environments,  with  buildings 
and  other  obstacles,  this  will  not  be  the  case.  Eor  this  situation  it  would  be  necessary 
to  incorporate  prior  knowledge  of  obstacles  (their  placement  and  materials  that 
determine  their  attenuation)  into  the  localisation  algorithms.  A  somewhat  similar 
but  different  problem  is  to  exploit  prior  knowledge  of  “forbidden”  zones  for  the  source 
location  (e.g.  secure  areas). 

3.  Diffusive  sources 
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Current  work  focused  on  point  sources  only.  In  some  situations  it  may  be  more 
realistic  to  use  a  model  of  spatial  and  temporal  concentration  distribution  of  the 
dispersed  substance  (chemical,  radiological)  from  a  diffusion  source  (e.g.  a  particle 
cloud  resulting  from  a  dirty  bomb  explosion).  The  goal  in  this  case  would  be  to 
track  the  centroid  and  the  spread  of  the  cloud  over  time.  In  the  case  of  a  chemical 
spill,  it  may  be  required  to  localise  the  source  of  diffusion. 

4.  Investigation  of  Wireless  sensor  network  technology  for  CBRN  data  fusion 

Recent  advances  in  wireless  sensor  networks  (WSNs)  technology  have  enabled  de¬ 
ployment  of  a  large  number  of  tiny  smart  sensor  nodes  for  monitoring  space  and 
objects.  The  unique  aspect  of  WSNs  is  that  they  integrate  wireless  communication, 
sensing  and  computation.  It  would  be  of  interest  to  investigate  the  use  of  tiny  CBRN 
sensors  (for  example  carried  by  soldiers)  and  their  fusion  for  monitoring  and  local¬ 
isation  of  CBRN  hazards  in  a  region  of  interest.  Various  aspects  of  this  emerging 
technology  would  need  to  be  examined,  such  as  the  fusion  architecture  (centralised 
vs  distributed),  network  bandwidth  and  latency  constraints,  protocols,  etc. 

5.  Optimisation  of  observer (s)  trajectory 

The  problem  is  to  guide  the  motion  of  (possibly  multiple)  observers  in  order  to  lo¬ 
calise  the  source  in  the  quickest  possible  manner  [17].  For  example,  suppose  that  we 
know  that  the  source  is  inside  a  certain  region,  and  we  perform  on-line  (recursive) 
source  position  estimation.  Then  intuitively  it  seems  that  the  observer(s)  should 
move  towards  the  currently  estimated  source  location,  in  order  to  get  more  informa¬ 
tive  measurements  and  localise  the  source  quicker.  The  observer  motion  should  be 
constrained  by  the  safety  of  observer(s).  Alternatively,  in  the  WSN  context,  we  may 
want  to  select  a  subset  of  sensors  to  transmit  the  measurements,  in  order  to  save  the 
batteries  (transmissions  require  energy  and  should  be  minimised  in  order  to  extend 
the  lifespan  of  a  network). 
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Appendix  A  Derivation  of  the  analytical  solution 


Let  us  consider  measurements  {zk,k=  1,2, 3, 4}  collected  at  four  arbitrary  points  {xk,yk)-  If 
we  disregard  noise  and  background  radiation,  these  measurements  can  be  expressed  as: 


2fc(x)=//r^,  A:  =1,2, 3, 4 


(Al) 


where 

Tk  =  \J {xk  -  xo)"^  +  ivk  -  UoY 

is  the  distance  from  the  source  {xo,yo)  to  the  fcth  measurement  point  {xk,yk)- 
These  four  equations  can  be  rewritten  as  follows: 

{xi  -  xof  +  {yi  -  yof  =  I j zx 

(x2  -  xof  +  (2/2  -  yof  =  Ilz2 

{xz  -  xof  +  (2/3  -  yof  =  I/Z3 

{X4  -  Xo)"^  +  (2/4  -  2/0)^  =  IjzA- 

By  subtracting  Eqn  (A4)  from  Eqn(A3)  we  get: 

{xl  -  xf)  -  2{xi  -  X2)xo  +  {yl  -  yi)  -  2(2/1  -  222)2/0  =  —  -  —  • 

Zl  Z2 


(A2) 


(A3) 

(A4) 

(A5) 

(A6) 


(A7) 


Similarly,  by  subtracting  Eqn  (A5)  and  Eqn  (A6),  respectively,  from  Eqn  (A3)  we  obtain: 


{xj  -  xl)  -  2{xi  -  xz)xo  +  {yl  -  yl)  -  2(2/1  -  2/3)220  =  —  -  — .  (A8) 

Zl  Zz 

and 

{xj  -  xl)  -  2(a:i  -  X4)xo  +  {yl  -  yl)  -  2(2/1  -  2/4)220  = - •  (A9) 

Zl  Z4 

By  dividing  Eqn  (A7)  by  Eqn  (A8),  we  get: 


{xl  -  xj)  -  2{xi  -  X2)xo  +  (2/1  -  yj)  -  2(2/1  -  222)2/0  ^  'k~ 
{xl  -  xf)  -  2(a:i  -  a:3)a:o  +  {yl  -  2/1)  -  2(2/1  -  2/3)220 

Similarly,  the  division  of  Eqn  (A7)  by  Eqn  (A9)  gives: 


(AlO) 


{xl  -  xj)  -  2(a:i  -  a:2)a:o  +  (2/1  -  yj)  -  2(2/1  -  2/2)220  ^  A  ~  7^ 
{xl  -  xl)  -  2{xi  -  a:4)a;o  +  {yl  -  yl)  -  2(2/1  -  2/4)220  ir  -  jr’ 

■^4 

Next,  we  define  the  following  variables: 

Ki  = 


{xl-xi)  +  {yl-yi) 
{xl-xl)  +  {yl-yl) 

{xl-xl)  +  {yl-yl). 


K2  = 

Li  = 
L2  = 
^3  = 


(All) 


37 


DSTO-TR-1988 


By  plugging  in  these  factors  in  Eqn  (AlO)  and  Eqn  (All)  and  rearranging,  we  obtain: 

aiXQ  +  biyo  +  ci  =  0  (A12) 

a2xo  +  b2yo  +  C2  =  0  (A13) 

where 


oi  =  2Ki{xi  —  X3)  —  2{xi  —  X2) 

bi  =  2Ki{yi  -  ys)  -  2{yi  -  y2) 

Cl  =  Li  —  K1L2 

02  =  2K2{xi  —  X4)  —  2{xx  —  X2) 

62  =  2K2{yi  -  yi)  -  2{yi  -  y2) 

C2  =  Li  —  K2L3. 

By  solving  Eqn  (A12)  and  Eqn  (A13)  simultaneously,  we  find  xq  and  yo  as: 

_  biC2  -  Cl 62 
0162  -  ^102 

_  O1C2  -  Cl  02 

61O2  -  0162  ’ 

By  inserting  the  values  of  Xq  and  j/q  in  Eqn  (A3),  we  compute  the  source  activity  I  as: 

I  =  zi  [(a;i  -  xof  +  {yi  -  yof]  ■ 


(A14) 

(A15) 

(A16) 
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